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I. INTRODUCTION 

Electronic correlations in materials have been one 
of the central topics of the condensed matter physics 
throughout its history encompassing topics such as 
high-temperature superconductivity, colossal magneto- 
resistance or heavy fermion physics. The introduction of 
dynamical mean-field theory (DMFT) in early 1990's^— 
marked a big step forward in the theory of correlated elec- 
trons in providing an approximate, but non-perturbative, 
computational method with several exact limits. More- 
over, numerical DMFT calculations are feasible also for 
multi-band Hamiltonians necessary for the description of 
real materials. 

In the past twenty years DMFT was applied, at first, 
to models^ and, later, to real materials^. Naturally, 
most of the studies focused on single-particle quantities, 
featured explicitly in the DMFT equations, which can be 
compared to photocmission spectra and which can pro- 
vide information about the metal-insulator transitions. 
Also local two-particle correlation functions, static as 
well as dynamical, can be computed with little additional 
effort in most DMFT implementations, providing infor- 
mation about the local response to applied fields. Com- 
putation of the non-local response is a more difficult ven- 
ture, however, with great potential gains. With the dy- 
namical susceptibilities available it is possible to compare 
to the experimental data from inelastic neutron, x-ray or 
electron loss spectroscopies. Perhaps even more inter- 
esting possibility opens with the static susceptibilities. 
Monitoring their divergencies as a function of tempera- 
ture and the reciprocal lattice vector allows investigation 
of the second order phase transitions and an unbiased 
determination of the order parameters. So far only cal- 
culations for simple models, such as the single-band Hub- 
bard model£i£,the periodic Anderson model^— , or Hol- 
stein modelii, have been reported. Similar calculations 
were also performed with cluster— or diagrammatic^ ex- 
tcntions of DMFT. 

In this paper, we present a scheme for computation of 
the static two-particle response functions in the multi- 
band Hubbard model within DMFT. The key develop- 
ment is splitting the particle-hole irreducible vertex into 
low frequency (LF) and high frequency (HF) parts, and 



expressing the HF part in terms of the local dynamical 
susceptibilities. This allows reformulation of the Bethe- 
Salpctcr equation in terms of the LF quantities plus cor- 
rections, thus reducing the numerical cost of the calcu- 
lations and improving their stability. The procedure is 
applied to the single-band Hubbard model at and away 
from the half filling and a two-band bilayer model. In all 
the studied cases we find a rapid convergence of the ver- 
tex function towards its asymptotic form, which leaves 
the size of the LF problem rather small and manageable 
also for multi-band systems. 

The paper is organized as follows. After a general in- 
troduction, the two-particle formalism is reviewed in sec- 
tion II, followed by the discussion of the numerical im- 
plementation. In section III applications to simple model 
systems are reported. Discussion of the asymptotic be- 
havior of the irreducible vertex and the blocked form of 
the irreducible Bethe-Saltpctcr equation is left to Appen- 
dices A and B. 



II. THEORY 

Our starting point is the multi-band Hubbard Hamil- 
tonian 

H = *R.-R' C R.i C R'j + 4 ^U,fc' c Ri c R.i c Rj c Rfc' 

<RR') R 

(1) 

where t l ^ i _- R , is the hopping amplitude between orbital 

j on site R' and orbital i on site R, (cr'j) are 
the corresponding creation (annihilation) operators, and 
Uij.ki is the anti-symmctrized local Coulomb interaction. 
Throughout the text we do not distinguish between spin 
and orbital degrees of freedom. Summation over repeated 
orbital indices is assumed. 

The developments and calculations reported here fall 
within the framework of dynamical mean-field theory. 
Detailed discussion of the formalism and basic applica- 
tions of DMFT can be found in Ref. 0. The main feature 
of DMFT is that the irreducible vertex functions (see be- 
low for details) are build only from local propagators and 
thus can be obtained from an effective quantum impurity 
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problem. Evaluation of the two-particle response func- 
tions can be viewed as a postprocessing of the solution to 
the DMFT equations, which self-consistently determine 
the fermionic bath for the effective impurity problem. 

A. Linear response formalism 

We review briefly the DMFT linear response formal- 
ism. For details the reader is referred to Ref. [j| We are 
interested in the response of a system controlled by fl) 
to static fields which couple to the spin, charge or a more 
general orbital density described by the susceptibility 

X«,«(q) = r^^e^ R ((T4.(r)c Rt (r)ct fc (0)c oi (0)) 

J ° R 

(2) 

Here q is a vector from the first Brillouin zone and the 
site indices were dropped in the local averages. In the 
DMFT approximation the susceptibility (0) can be ob- 
tained from 

Xy,jw(q) =T^2xijM c i>< J ™> UJ n ), (3) 

where Xij,fc/(q; w m , oj n ) is the solution of coupled integral 
equations^ 

x l j,fcz(q;wi,w 2 ) = x°j,ki (q; <•>!., w 2 ) 
+ t 53 

(4) 

Xij,*i(wi,W2) = x%,ki{ux,u 2 ) 

+ T ^Z X°j,mn( w l7W3)r m „^ 9 (w3,W 4 )Xp I? ,A ; ;( w 4,W2)- 

(5) 

The local and corresponding 'lattice' quantities are dis- 
tinguished by the presence of parameter q. The quanti- 
ties is the equations arc functions of the discrete fermionic 
Matsubara frequencies, which arc at temperature T given 
by uj n = (2n — 1)tvT (n is an integer). The equations 
have the form of the irreducible Bethe-Salpeter (IBS) 
equation, with r mn>pq (uj3, 0J4,) being the particle-hole ir- 
reducible vertex at zero energy transfer, in the diagram- 
matic expansion of which the run pair of external ver- 
tices cannot be separated from pq pair by cutting two 
single-particle propagators. In the DMFT approximation 
only the local diagrams contribute to T mniPq (u>3,ui4)^, 
which allows simplifying the IBS equation for the full 
two-particle correlation function to equation in terms 
of the reduced quantity Xij,jy(q;k>i, w 2 ) and k-integrated 
particle-hole bubble 

X^-,fc;(q; w i> w 2) = -^jp-J2G lk (k]UJ 1 )Gi j (k + q;uj 1 ). 

k 

(6) 



Gifc(k; uji) is the single-particle propagator obtained from 
the Dyson equation 

G ifc (k;w n ) = [iw n + n - h k - Jl(uj n )]]. k , (7) 

where fj, is the chemical potential, h k is the Fourier trans- 
form of the single-particle part of the Hamiltonian (TTJ) 
and £(w„) is the local self-energy. Equation ([5]) is the IBS 
equation for the impurity, relating the local two-particle 
correlation function 

rP rP rP rP 
Xij,fei(wi, W2) = T 2 / dn / dr 2 / dr 3 / dr 4 
Jo Jo Jo Jo 

e iWl (— 2 ) e - 2 ^-^((Tc i (r 1 )c)(r 2 )c,(r 3 )4(r4)) (8) 
-(Tc i (r 1 )ct(r 2 ))(Tc J (r 3 )4(r4))) 

to the local particle-hole bubble 

X°j,fei(wi,w 2 ) = -^^^Gikfcuj^Gijik';^). (9) 



B. Numerical implementation 

Although it is possible to eliminate the vertex function 
from equations (j4j and ([5]), we evaluate Tij ^1(011,102) ex- 
plicitly as an intermediate step. The computation pro- 
ceeds as follows. After converging the DMFT equa- 
tions we use the impurity self-energy Hij(ui n ) to get 
X°^fei(q; w i,w 2 ) and X°j,fei(^i)W 2 ) using expressions @, 
(0), and (J7J. The local two-particle correlation function 
(wi,o; 2 ) is obtained during the solution of the im- 
purity problem. Next, equation ([5]) is inverted to obtain 
ry,fci(wi,a; 2 ). Finally, equation ((4]) is solved and the sus- 
ceptibility Xij.fcz(q) obtained after summation over the 
Matsubara frequencies ([3]). 

The above program faces two numerical challenges. 
First, the inputs to the IBS equations, Xy,fci(wi,Cc> 2 ) in 
particular, are known accurately only for a limited range 
of frequencies. Second, Xii.fcK^l:^), X%,klfai u l> 
Xij,fci(q; wi,w 2 ), and ^.^(q; Wi, w 2 ) decay as 1/w 2 , 
which makes straightforward inversion of ((4]) and (|5|) nu- 
merically unstable. We have been using the Hirsch-Fye 
quantum Monte- Carlo impurity solver—, but these issues 
are general and apply to other methods as well. Our com- 
putational procedure relies on splitting the problem into 
a low and a high frequency part, which arc treated in 
numerically different ways. 



1. Asymptotic behavior of the vertex function 

Separation of low and high frequency parts of a given 
problem is common in theoretical physics. In the DMFT 
practice it is widely used in calculation of the self-energy, 
where the numerical solution of the Dyson equation at 
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low-frcqucncics is augmented by a high-frequency asymp- 
totic expansion obtained from the moments of the spec- 
tral function^. In Appendix A we start by considering 
the HF asymptotic expansion for the self-energy, which 
serves a precursor for the HF expansion of the vertex 
function IY^^i, wa)- 

Unlike the self-energy, where a separate equation ex- 
ists for each frequency and so the LF and HF parts are 
decoupled, in case of equations dU and §5§ we are dealing 
with matrices in the Matsubara frequencies and have to 
solve coupled equations for the HF and LF blocks. The 
key ingredient of our procedure is replacing the vertex 
function outside the LF block by its asymptotic form 
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w 2) = UijM + f™ ; hXmV w ( u l - W2)U v j A l 

for > uj c V |w 2 | > cj c , 

(10) 

where x ph { v ) and X pp (f) are the local dynamical suscep- 
tibilities in the particle-hole and particle-particle channel, 
respectively, as functions of bosonic Matsubara frequency 
v, and Lu c is the cut-off frequency defining the LF block. 



dr exp(i^) ( (Tcj (r)c t (r)c{ (0)c, (0) ) 



(11) 



X^«M= / *-«xp(ii^)(Tc i (r) Cj .(T)4(0)4(0)> (12) 



FIG. 1: The Neel temperature of the 2D Hubbard model 
in DMFT approximation (squares), shown together with the 
mean-field value of the Heisenberg model (red, large U) and 
the RPA results obtained with the non-interacting (dotted 
blue line) and dressed (triangles) propagators. The inset 
shows the decrease of the Neel temperature upon doping away 
from half-filling at U — 1. 



where each element x is itself a matrix in pairs of orbital 
indices 
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11 / Xll.ll Xll,21 
21 X21,ll X21,21 



X = 



12 



Xl2,ll Xl2,21 
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In this representation the LF block is located in the upper 
left corner of the matrix and the uj dependent part of 



Derivation of expression ©is given in Appendix A. p«, ^ m proportional to X ph and X pp has a band 



ij,kl 

constant, U^m, plus two ridges along the main and the 
minor diagonal, a structure observed in earlier numerical 
studies^. The cross section of these ridges is given by 
the local dynamical susceptibility, typically, with a fast 
decay away from v — 0. The asymptotic form (|10[) is 
used in the numerical treatment of both equations (|4I5[) 
as described below. 



2. Blocked IBS equations 

The relations (|4l5p can be viewed as equations between 
matrices indexed by the Matsubara frequencies and pairs 
of the orbital indices. We use the following arrangement 
of the matrices 



x(+^i,-^2) •• 

X(-W2,-Wl) X(-W 2 ,+Wi) X{-U2,-U2) ■■ 



diagonal form. 

In the following we split the set of Matsubara frequen- 
cies into the LF block |w n | < u> c , denoted with a block 
index 0, and HF block \ui n \ > w c , denoted with a block 
index 1 . Expansion of a matrix equation of the form (|4I5[) 
into the LF and HF blocks is given in Appendix B. The 
impurity IBS equation (0 is used to compute the LF 
part of the vertex function, the T 00 block. Wc assume 
to know the LF (00) part of Xy.fci^ii^a) © from the 
solution of the impurity problem and to have a complete 
information about %° fc; (wi,cj 2 ) ©■ In practice we use 
a numerical representation of x% fei( w i! ^2) up to a HF 
cut-off f2 c (typically much larger that the LF cut-off ui c ) 
and the analytic — -\ tail up to infinity for elements the 



X 



ij-ij 



elements, while the remaining elements involving 
off-diagonal G%j §2§ arc neglected above fl c . The LF part 
of the vertex T 00 is obtained from (|B4| , which amounts 
to inverting equation ([5]) restricted to the LF block and 
subtracting a correction. For x vh { v ) and x pp (^) sharply 
peaked around v = the correction reduces to a constant 
shift everywhere except the vicinity of w c . 

Once we have computed T 00 wc proceed to solve equa- 



FIG. 2: The real part of the transverse spin vertex 
ri2i2 (k> m , u) n ) as a function of the indices m, n for U = 0.4 
(left) and 1 (right) obtained at T = 1/15. 
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FIG. 3: The real part of the transverse spin vertex 
Ti2i2 (^m, <^n) for fixed u) m in the asymptotic range for U = 
0.4 (left) and 1 (right) at T = 1/15. Full diamonds mark the 
numerical data from LF block, the crosses denote the HF tail 
obtained from expression 1101 



tion ([I]) for Xij,kl(fa<*>i>i*>2) assuming the full knowledge 
of X° ifc/ (q; wi,<j 2 ) and ry^jfwijWa). We use the same 
partitioning into the LF and HF blocks and take ad- 
vantage of the fact that only w-summed susceptibility 
Xij,fcz(q) is °f interest. This allows partial ^-summations 
in the HF block to be preformed at an earlier stage of the 
calculation and to replace matrix-matrix operations with 
matrix-vector ones, see Appendix B. Further speed-up 
comes from treating the cj-dependent (|B11[) and the con- 
stant (|B12[) parts of T 11 in two subsequent steps, which 
reduces the most demanding part of the HF problem to 
a repeated band-matrix-vector multiplication. 



III. NUMERICAL RESULTS 

A. 2D Hubbard model 

In the following we present our results for 2D Hubbard 
model on a square lattice with the nearest neighbor hop- 
ping. In particular, we compute the Neel temperature as 
a function of U and doping. This is an academic example 
since it is well know that no magnetic order is possible in 
this model at finite temperature and the magnetic oder 
studied here is an artefact of the DMFT approximation. 
Nevertheless, we find it a useful benchmark since i) a 
comparison to other studied is possible and ii) the model 
exhibits a crossover from a nesting-driven instability at 




FIG. 4: The longitudinal spin susceptibility as a function of q 
at various temperatures for U = 0.4 (left) and U = 1 (right). 



small U to the Heisenberg model with local moments at 
large U. 

First, we report the results obtained at half filling. Af- 
ter converging the DMFT equations, a finer QMC run of 
5 x 10 6 to 10 7 sweeps is performed to obtain the local two- 
particle correlation function Xij,kl{wi, 0J2) and local sus- 
ceptibilities x P ijki^ u )^ X P ijki^ v )- It i s possible to separate 
the charge and the spin-longitudinal channels from the 
beginning, but we perform the calculation in the upup, 
updn, dnup, dndn basis, as would be the case for general 
orbital indices, and separate the two channels only in the 
final result. Expression ^ is evaluated on a 20 x 20 grid 
of q-points (66 irreducible points) using a fine k-mesh of 
551 x 551 points, which is more than enough to ensure 
that the uncertainties associated with k-points sampling 
as much smaller that other sources of error. 

The IBS equations were solved for several LF cut-offs 
u> c leaving between 5 to 20 lowest Matsubara frequencies 
for temperatures above 1/20, and 15 to 30 Matsubara fre- 
quencies for lower temperatures. In all studies cases, the 
relative variation of the susceptibility with the LF cut-off 
was less than 1%. We have also tested the sensitivity to 
the HF cut-off il c . When neglecting the contributions 
from above £l c completely, a slow convergence of the sus- 
ceptibilities close to the AFM instability was observed. 
The error essentially followed the l/f2 c dependence which 
may be expected when truncating a sum over se- 
ries. The dominant contribution to the error came from 
the r 01 X 11 r 10 term in equation (|B6|) . Introduction of 
an analytic summation of the l/w„ tails to infinity lead 
to the converged results already for lowest HF cut-offs 
corresponding to 80-100 lowest Matsubara frequencies. 

To check the consistency of our calculations we have 
solved equation Q also for the local bubble Xij u ( Wl ' W2 ) > 
i.e. went from Xij,kl( u it ^2) to Ti^ui^u^i) and back, 
and compared the result with the local susceptibility ob- 
tained directly from the QMC simulation. For all re- 
ported parameters we have found a relative deviations of 
less than 2% close to the AFM instability and less than 
0.5% deeper in the paramagnetic phase. As a side remark 
we mention that in the insulating phase for U > 1 it is 
important that Xij ki&i'^z) entering equations Q and 
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(([5]) are exactly the same (e.g. obtained with the same 
k-point sampling) since in this parameter range the re- 
sults are extremely sensitive to differences between x°'s 
entering equation ^ and ([5]). 

Due to the SU(2) symmetry of the model the transverse 
and the longitudinal spin susceptibilities are identical in 
principle. This symmetry is not explicitly enforced in the 
QMC calculation, where the identity is obeyed only ap- 
proximately within the stochastic error of the simulation. 
For U < 1 we have found relative deviations between the 
transverse and longitudinal spin susceptibilities of 0.2- 
0.5%. For larger f/'s (especially at lower temperatures) 
the transverse susceptibility rapidly becomes useless due 
to increasingly large errorbars. This problem is specific to 
the Hirsch-Fye QMC solver and the way the transverse 
susceptibility is measured in the simulation. We have 
used the longitudinal susceptibility to determine Tm in 
this parameter range. 

In Fig. Q] we show the Neel temperature as a func- 
tion of U. On the large U side Tjv approaches the static 
mean field value of the corresponding Hcisenberg model 
reflecting the insulator with well defined local moments. 
On the small U side the magnetic order arises from the 
Fermi surface nesting. For comparison, we show also T/v 
for the RPA condition for the instability, x°(T/v) = 1/U, 
obtained by replacing the vertex with a constant and tak- 
ing x°(T i ) of the non- interacting electron gas or x°(T) 
obtained with the dressed propagators. 

In Fig. [2] we show the LF block of the real part of 
the transverse spin vertex r 12i i2 for large and small U 
as a function of the fcrmionic Matsubara frequencies, 
which exhibits the two-ridge structure and a rapid con- 
vergence towards its asymptotic form. At small U the 
ridges arc only weakly T-dcpcndent and their absolute 
value is small compared to the constant part of the ver- 
tex. At large Us the ridges grow as 1/T and thus be- 
come dominant at low enough temperatures. In Fig. [3] 
we show the same vertices at fixed u> m (still in the LF 
block) together with the corresponding asymptotic tails 
rjSj 12 . We find a good match between the two depen- 
dencies. The remaining mismatch is mainly due to the 
time discretization error inherent to the Hirsch-Fye im- 
purity solver, which shows up in a spurious curvature 
of the T{uj m ,uj n ) for larger ws (not really visible in the 
present plots). Using the semi- analytic asymptotic tail 
T°° can efficiently remove this spurious curvature. The 
discretization error arises from performing Greens func- 
tion convolutions on a discrete time grid and should not 
be confused with the so called Trotter error, which is 
another consequence of discretizing the imaginary time 
and which cannot be removed be restricting oneself to 
the low frequencies. Impurity solvers using continuous 
time approaches are not expected to have this problem. 

Finally, the spin susceptibility throughout the Bril- 
louin zone is shown in Fig. [4] The plots clearly show 
that the magnetic instability takes place at (tt, tt) as ex- 
pected. The behavior of the vertices reveals that the ori- 
gin of the instability is quite different for small and large 
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FIG. 5: The detail of the q-dependent spin susceptibility in 
the vicinity of [tt, tt] at T = 1/60 (larger values) and T = 1/30 
for the doping of 0.18. 

Us. While for large Us the main T-dependence in the 
problem is in the vertex, related to 1/T behavior of the 
local susceptibility, for small Ux the main T-dependence 
is in the 'bubble', namely the of thermal smearing of the 
nesting related peak. 



B. Effect of doping 

We use doping to break the particle-hole symmetry in 
the above model. The calculations were performed in 
the crossover regime between the metallic and insulating 
phases at U = 1. We have found that magnetic insta- 
bility survives doping up to about 0.2 electron (or hole) 
per site. For dopings less or equal to 0.14 the instabil- 
ity appears at [7T,7t]. For dopings of 0.16 and 0.18 we 
have noticed that the maximum of Xs (q) moved slightly 
away from [tt, tt] as the temperature was lowered. Simi- 
lar behavior was observed also in a model with the next- 
nearest-neighbor (diagonal) hopping^. The Neel tem- 
perature as a function of doping is shown in the inset of 
Fig. Q] In Fig. [5] we show the detail of the spin suscepti- 
bility in the vicinity of [tt, tt] for the highest studied dop- 
ing of 0.18. While at high temperatures the susceptibility 
exhibits a flat maximum at [tt, 7r], at lower temperatures 
the maximum moves to [1.047r,7r] eventually giving rise 
to an incommensurate order. 



C. Bilayer model 

As a last example we study a bilayer Hubbard models. 
We choose this model because: i) it has two bands and 
the lattice bubble Xij,kl(m W\,ui2) has a more complicated 
orbital structure than for the single-band model ii) it ex- 
hibits a non-trivial temperature dependence of the uni- 
form spin susceptibility characterized by an exponential 
decrease at low temperatures. In our previous work on 
a more general version of this model (model of covalent 
insulator )2£ we have used a different numerical approach 
to compute the uniform susceptibility without an explicit 
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IV. CONCLUSIONS 
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FIG. 6: The uniform spin (circles) and charge (squares) sus- 
ceptibility of the bilayer model as a function of temperature. 
The diamonds mark the local spin susceptibility, which clearly 
shows the evolution of local moment above T = 0.1. 



An efficient numerical procedure for calculation of the 
static two-particle response functions within the DMFT 
formalism was presented. The main development is the 
identification of the asymptotic form of the two-particle 
irreducible vertex function and implementation of a block 
algorithm for solution of the irreducible Bethe-Salpctcr 
equations. We have studied several test cases includ- 
ing the single-band Hubbard model at and away from 
half-filling, and a bilayer Hubbard model at half-filling 
using Hirsch-Fye quantum Monte-Carlo impurity solver. 
In all cases we have observed a rapid convergence of 
the vertex function towards its asymptotic form and a 
very good numerical stability of the computations. This 
means that the numerical effort can be efficiently re- 
duced by choosing a relatively small low-frequency cut- 
off, which makes calculations for multi-band systems 
feasible. Besides the susceptibility calculations the de- 
scribed procedure may be useful in diagrammatic exten- 
sions of DMFT, which employ the IBS formalism and 
two-particle vertice s 17 ' 21 " — . 



use of the vertex function and its HF tail, an approach 
which could not be used for more general models and 
q-dependent susceptibilities. Obtaining the exponential 
decay at low temperatures proved then numerically chal- 
lenging. Therefore we find this model to provide a good 
test of the performance of the present method. 

The model is obtained by taking two Hubbard sheets 
from the previous section (with opposite signs of the 
in-plane hopping) and adding an inter-plane hopping of 
0.2. The different signs of the in-plane hoppings on the 
two sheets ensures that any non-zero inter-plane hopping 
opens a hybridization band gap. The inter-plane hopping 
was chosen such that the interesting temperature range 
can be easily covered with the QMC simulations. As 
in the other cases studied so far we were not primarily 
interested in the physics of the model, but in the per- 
formance of the computational method. Therefore we 
omit that fact that the system actually exhibits an AFM 
instability and focus on the uniform susceptibility only. 
While the AFM instability can be removed by introduc- 
ing an in-plane frustration, the T-dependence of the uni- 
form susceptibility is a generic feature independent of the 
bandstructure details^. 

In Fig. [6] we show the local and uniform spin suscep- 
tibilities together with the uniform charge susceptibil- 
ity, which is inversely proportional to the compressibility 
of the electron liquid. The calculated spin susceptibili- 
ties closely resemble those of Ref. HtJ The decay of the 
uniform susceptibilities at low temperatures reflects an 
evolution towards a band-insulator-like state, character- 
ized by a finite excitation gap for both charge and spin 
excitations (for more information see Ref. Il9i ). In the 
present calculations we have not encountered any nu- 
merical problems in reproducing the low values of the 
susceptibility at low T. 
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Appendix A: High frequency expansion of the 
vertex function 

Following Baym and KadanoflS^ and Baym2£ the self- 
energy S can be viewed as a functional of the dressed 
Green function G and the vertex function T can be ob- 
tained as a variation ■ The self-energy itself can 
be obtained as a variation of the Baym-Kadanoff func- 
tional E = • Diagrammatically this means that E is 
obtained from $ by cutting a single Green function line 

and r = S gQg^ is obtained by cutting two Green function 
lines. 

We are interested in the high frequency behavior of 
E(w) and r(wi,W2)- The expansion of E(w) to the or- 
der ^ is well known from the moments of the spectral 
function-^. Diagrammatically the zero's order term of E 
is the Hartrec-Fock contribution, in which the external 
frequency does not enter any of the internal lines (Fig. 
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(a) (b) 

FIG. 7: The form of the diagrams for the self-energy S(^), 
which contribute to the zero's (a) and first (b) order in i. 




(a) (b) (c) 

FIG. 8: The diagrams for the irreducible vertex, which remain 
non-zero in the limit of large uj\ and 0J2 ■ 

[7Ji). The first order contribution, 1/uj, comes from dia- 
grams with external vertices connected by single Green 
function line (Fig. Wp)- Their contribution has the form 

T^UxMUG^-u) (Al) 

with xW) being the dynamical susceptibility, where the 
plus sign in front of u applies for particle-hole and the 
negative sign for particle-particle susceptibility. Expand- 
ing G in (|A1[) in the powers of j-^p. = ^ + ZiU^lT) an< ^ 
taking the leading — term we get the corresponding con- 
tribution to the self-energy. The remaining free sum- 
mation over v leads to an equal-time correlator of the 
(c^c c'c ) type. 

Our main interest here lies in the high frequency be- 
havior of T(lui,uj2). In particular, we are interested in 
the zero's order contribution, which remains finite in the 
limit or large u>i and/or u)2- It is not possible to de- 
velop the high frequency expansion by keeping one of 
the frequencies fixed since the limit u>i , uj-2 — > 00 is not 
uniform, i.e. it depends on the way in which infinity is 
approached. By inspecting the diagrams for T(u>i,u> 2 ) 
we have identified two non- vanishing contributions. The 
first one, obtained from the variation of the zero's order 
term in the self-energy (Fig. [7^,) is simply the bare vertex 
Fig. [Sh*. The other non-vanishing terms are obtained by 
cutting the Green function line in Fig. [TJd, which leads 
to Fig. |SJd and Fig. [5J:. The corresponding expression 
for r°° reads The expression r°° reads 

+ ^Uim^iXmn.pqi^ 1 + ^Wp^kq ■ 

(A2) 



(a) (b) (c) 

FIG. 9: The diagrams in the Baym-Kadanoff functional $[G], 
which contribute to F°° . 

The factor 1/4 accounts for the exchanges of the external 
points of the two vertices. Note, that the second order 
diagram (obtained by replacing the susceptibility in Fig. 
[7b with a simple bubble) contributes to both Fig. [5b and 
Fig. [SJ;. As a function of uj\ and uj 2 the high frequency 
part of the vertex function r°° has a structure of a con- 
stant plus two ridges along the main and minor diagonals, 
cross section of which is determined by the particle-hole 
and particle-particle dynamical susceptibilities. 

In terms of the Baym-Kadanoff functional, the con- 
stant term in T°° derives from the first order, Hartree- 
Fock, diagram (Fig. E^), and the diagrams containing 
at least one particle-hole or particle-particle bubble (Fig. 
[5Jd,c) by cutting the Green function lines corresponding 
to the same bubble. 



Appendix B: Blocked IBS equations 

In the following we present the IBS equation written in 
terms of small-w and large-w sectors denoted by indices 
and 1, respectively. The IBS equation has the form 

A = B + BTA, (Bl) 

where A, B and T arc matrices indexed by the discrete 
Matsubara frequencies. Matrix B is diagonal in the Mat- 
subara frequencies. 

In the first step we assume that A 00 , B 00 , B 11 , T n , 
r , and r 10 are known and we want to compute T . 
From 

A 00 = b 00 +B 00 T 00 A 00 + B m T 01 A w (B2) 
A w =B n r io A oo + s ii r ii^io (B3) 

we obtain 

A oo = B oo + ^oo^oo + r oi x ii r io> )A oO ) (B4) 

where X 11 full-fils the equation 

X 11 = B 11 + B 11 T 11 X 11 . (B5) 

Thus r 00 can be obtained by inverting IBS equation trun- 
cated to the 00 block and adding a correction term. 

In the second step we assume full knowledge of T and B 
and want to compute A. Importantly, we are interested 
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only in the sum over the matrix elements (Matsubara 
frequencies) of A. To this end we use a set of block 
equations 

A ao = B ao + B oo( r oo + pOi^iipio^oo (Bg) 

A 10 = x n r 10 ^ 00 (B7) 
A 01 = A 00 r 01 X n (B8) 
A 11 = X 11 + X^A^T^X 11 . (B9) 

The advantage of the blocked scheme is that equation 
(|B5|) docs not have to be inverted in practice. This is 
thanks to three factors: i) only sum of elements of X 11 
(or a weighted sum in case of X^T 10 ) is needed, ii) B 11 
is small due to 1/w^ decay, iii) T 11 is a constant plus a 
band matrix. We describe in detail the computation of 

x. m = J2 x a- ( B1 °) 

i,3 



The computation for X n T w is analogous. First, we solve 
iteratively equation (|B5[) taking only the band part of T 11 

Xi. = B u + B a r^ nd X k ., (Bll) 

where the block indices 11 were dropped for simplicity. 
After X^ is computed the summation over the remaining 
index is performed to obtain X, m . Next, we should solve 
equation (|B5|) taking the constant part of T 11 and X 
in place of B 11 . This is not necessary, since there is a 
simple, well known, relationship between X„ and X %% 

X.. = (/ - l..r const ) _1 l... (B12) 

Note, that quantities in this equation should be viewed 
as non-commuting matrices in the orbital indices. 
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